package org.evolver.disser;

import org.apfloat.Apfloat;

public class RootFinder {
	private final int		MAX_ITERATIONS    = 20;
	private final int		DEFAULT_PRECISION = 75;
	private final Apfloat	APFLOAT_TWO       = new Apfloat(2.0, DEFAULT_PRECISION);
	private final Apfloat	APFLOAT_THREE     = new Apfloat(3.0, DEFAULT_PRECISION);

	/**
	 * 
	 * @author evolver
	 *
	 */
	private class Interval {
		public Interval() {
			from = null;
			to   = null;
		}

		public Interval(Apfloat fromIn, Apfloat toIn) {
			from = fromIn;
			to   = toIn;
		}

		public Apfloat	from;
		public Apfloat	to;
	}
	
	/**
	 * 
	 * @param func
	 * @param from
	 * @param to
	 * @return
	 */
	public Apfloat FindRoot(Function func, Apfloat from, Apfloat to) {
		Interval	adjustedInterval = AdjustInterval(func, from, to);
		if (adjustedInterval == null) {
			// This is a fatal error. Bail out.
			return null;
		}
		
	    Apfloat xLeft     = adjustedInterval.from;
	    Apfloat yLeft     = func.GetValue(xLeft);
	    int     signLeft  = yLeft.signum();

	    Apfloat xRight    = adjustedInterval.to;
	    Apfloat yRight    = func.GetValue(xRight);
	    int     signRight = yRight.signum();

		//  A very special case when both boundaries are roots and yet
		//  they have different values. This should never happen but
		//  for genericness we have to support that.
	    if ((signLeft == signRight) && (signLeft == 0) && !(xLeft == xRight)) {
	        return null;
	    }

	    //  Special case when left boundary is a root
	    if (signLeft == 0) {
	        return xLeft;
	    }

	    //  Special case when right boundary is a root
	    if (signRight == 0) {
	        return xRight;
	    }

	    //  See if root can be found. 
	    if (signLeft == signRight) {
	        return null;
	    }

	    Apfloat     xMiddle;
	    Apfloat     yMiddle;
	    int         signMiddle;
	    for (int i = 0; i < MAX_ITERATIONS; ++i) {
	        xMiddle    = xLeft.add(xRight).divide(APFLOAT_TWO);
	        yMiddle    = func.GetValue(xMiddle);
	        signMiddle = yMiddle.signum();
	        if (signMiddle == signLeft) {
	            xLeft = xMiddle;
	            continue;
	        }

	        if (signMiddle == signRight) {
	            xRight = xMiddle;
	            continue;
	        }

	        xLeft  = xMiddle;
	        xRight = xMiddle;
	        break;
	    }

		return xLeft.add(xRight).divide(APFLOAT_TWO);
	}

	/**
	 * 
	 * @param func
	 * @param from
	 * @param to
	 * @return
	 */
	public Interval AdjustInterval(Function func, Apfloat from, Apfloat to) {
		Interval	ret = new Interval(from, to);

		//  Let's set the initial offset at 1/3 of the original interval length
		//  It has to be less than 1/2 of it, though.
	    Apfloat offset = to.subtract(from).divide(APFLOAT_THREE);

	    for (int i = 0; i < MAX_ITERATIONS; ++i) {
	        Apfloat xLeft    = from.add(offset);
	        Apfloat yLeft    = func.GetValue(xLeft);
	        int     signLeft = yLeft.signum();

	        Apfloat xRight    = to.subtract(offset);
	        Apfloat yRight    = func.GetValue(xRight);
	        int     signRight = yRight.signum();

	        // Let's make the offset smaller...
	        offset = offset.divide(APFLOAT_TWO);;

	        // Special case when left boundary is a root
	        if (signLeft == 0) {
	            ret.from = xLeft.subtract(offset);
	            ret.to   = xLeft.add(offset);
	            return ret;
	        }

	        // Special case when right boundary is a root
	        if (signRight == 0) {
	            ret.from = xRight.subtract(offset);
	            ret.to   = xRight.add(offset);
	            return ret;
	        }

	        // See if root can be found on the current interval
	        if (signLeft != signRight) {
	        	// Yes, it can!
	            ret.from = xLeft;
	            ret.to   = xRight;
	            return ret;
	        }
	    }

	    //  If we reached that far, we could not stretch the interval
	    //  big enough to enclose a root.
	    return null;
	}
}
